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1.  Introduction 


The  Army  Research  Laboratory  (ARL),  Atmospheric  Laser  Optics  Testbed  (ALOT)  is  a  unique 
experimental  facility  with  which  to  measure  optical  turbulence  ( Cn2 )  intensity  and  investigate  its 
effects  on  infrared  imaging  and  laser  optics  communications  (1,2).  Within  the  A  LOT,  a  near 
horizontal,  2.3  km  optical  path  extends  from  the  top  of  a  tall  water  tower  to  the  Intelligent  Optics 
Laboratory  (IOL)  rooftop  at  ARL.  However,  complex  microphysical  influences,  e.g.,  irregular 
wind  flow  patterns  around  the  IOL  and  the  water  tower  and  the  effects  from  wind  shears  and 
temperature  changes  across  the  top  of  nearby  forest  canopies,  may  affect  the  A  LOT  measured 
data  and  research  applications.  To  this  end,  computer  simulation  models  may  provide  some 
meaningful  results  even  though  all  the  pertinent  landscape  and/or  canopy  characterization  data 
along  the  optical  path  may  not  yet  be  known  or  available. 

In  a  previous  study,  Tunick  (3)  reported  on  the  characterization  of  complex  (e.g.,  urban)  signal 
propagation  environments  via  physics-based,  computational  fluid  dynamics  (CFD)  models. 
However,  CFD  program  codes  are  (as  a  rule)  quite  computationally  intensive  (4,5).  For 
example,  CFD  codes  can  require  1  to  8  hours  or  more  of  execution  time  on  multiprocessor  super¬ 
computers.  In  addition,  CFD  models  are  generally  cumbersome  to  modify  and  debug.  In 
contrast,  it  was  found  that  a  useful  mathematical  representation  of  several  key  microphysical 
parameters  inside  and  above  forests  could  be  obtained  by  means  of  a  rapid,  steady-state,  second- 
order  turbulence  closure  model,  with  an  embedded  radiative  transfer  and  energy  budget 
algorithm  to  predict  the  heat  source  (6, 7).  However,  development  of  this  later  computer  model 
made  it  possible  to  generate  realistic  profiles  for  the  wind  speed  and  temperature  within  and 
above  uniform  vegetative  canopies  only.  Hence,  we  desired  an  intermediary  model,  which 
would  be  computationally  efficient  and,  at  the  same  time,  complete  in  physics  to  characterize  the 
propagation  environment  around  different  building  geometries  embedded  in  the  model  grid  (as 
well  representation  of  canopy  drag  forces).  In  other  words,  we  are  interested  in  models  that  are 
still  reasonably  fast  but  have  enough  flexibility  to  apply  to  the  types  of  field  tests  envisioned  in 
future  works  (e.g.,  studies  related  to  free-space  laser  optics  communications  from  within  the 
ARL  A  LOT). 

In  this  paper,  we  present  a  finite-difference  computer  model  to  predict  the  three-dimensional 
microphysical  influences  on  signal  propagation  in  complex  areas,  e.g.,  around  single  and 
multiple  building  arrays  and  forests.  So  far,  the  model  accounts  for  advection,  the  pressure 
gradient,  and  drag  forces  due  to  vegetation  (e.g.,  open  fields  or  forests).  The  model  incorporates 
first-order  (eddy  diffusivity)  turbulence  closure  for  neutral  stability  (see  Mellor  and  Yamada  (S)). 
Note  that  mechanisms  to  account  for  heating,  cooling,  and  moisture  flux  are  not  yet  considered 
(in  order  to  maintain  computational  efficiency  and  flexibility  with  regard  to  modifications  and 
debugging).  Nevertheless,  our  research  strategy  includes  the  development  of  such  algorithms, 
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which  can  be  implemented  in  future  works  to  improve  existing  optical  turbulence  ( Cn2 ) 
modeling  and  analysis  around  the  ARL  A  LOT  Facility. 


2.  Model  description 


Following  the  numerical  methods  discussed  by  Meyers  (9)  as  well  as  those  found  in  several  key 
texts  on  fluid  mechanics  (10-14),  the  simplified  Navier  Stokes  equation  (in  non-conservative 
form)  for  the  current  model,  neglecting  Coriolis1  forces,  can  be  expressed  as 

^  +  u  •  Vn  =  J3v(y  ■ u)+kV2u+f  ( 1) 


where  t  is  the  independent  variable  time,  u  =  u,  v,  w  is  the  wind  velocity  vector,  ft  is  a  constant,  k 
is  the  eddy  viscosity  (assumed  constant),  and  /  is  the  drag  force  acceleration  due  to  vegetation. 
In  equation  1 ,  the  advection  term  is 


and 


(u  •  V  )m  = 
(w  •  V  )v  = 


r  d  d  d A 

u - 1-  v - Mv — 

dx  dy  dz 


d  d  d 
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dx  dy  dz 


(u  ■  V  |vv  ■ 


r  d  d  d  A 
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^  dx  dy  dz  j 
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(2b) 


(2c) 


The  divergence  portion  of  the  pressure  gradient  term  is 


.  -  ^  du  dv  dw 

A  =  V  -u=  —  +  —  +  — . 

dx  dy  dz 


(3) 


and  the  Laplacian  portions  of  the  eddy  diffusion  term  are 


v/2 

V  u  - 


o  u 


yi 
o  u 


dx  dy i 


d2u 

5? 


(4a) 


^  An  earlier  survey  of  the  literature  (IS)  noted  that  other  than  a  few  authors,  e.g.,  Shinn  (16),  Holland  (1 7),  and  Wilson  and 
Flesch  (18),  most  exclude  the  effect  of  the  Coriolis  force  as  having  negligible  effect  on  the  scales  of  motion  considered,  i.e.,  wind 
flow  through  the  forest  canopy  layer. 
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(4b) 
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(4c) 


Finally,  the  vegetation  drag  force  accelerations  are 

fx  —  CjA\u\u ,  (5a) 

fy^CdA\u\v,  (5b) 

and 

fz=cdA  \"\w>  (5c) 

where  Q  is  the  drag  coefficient  (0. 1-0.2),  A  is  the  leaf  area  density  (see  ref.  ( 6 )),  and  \u\  is  the 
magnitude  of  the  total  wind. 

Equation  1  is  solved  explicitly  forward  in  time,  with  the  advection  term  upwind  differenced  in 
space  and  the  diffusion  term  central  differenced  in  space.  As  an  example,  the  equation  for  the  u- 
component  of  the  wind  can  be  discretized  as 
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(6) 


The  boundary  conditions  applied  in  this  model  are  as  follows:  inflow  (typically  from  the  west)  is 
a  specified  value  (Dirichlet  boundary  condition),  such  as  a  constant  or  log-law  velocity  profile. 
Outflow  (typically  to  the  east)  is  constrained  to  have  a  constant  flux,  i.e.,  zero  gradients 
(Neumann  boundary  condition).  The  north  and  south  sides  are  also  constrained  to  have  a 
constant  flux.  The  vertical  velocity  at  the  bottom  of  the  model  is  a  specified  value  (e.g.,  w  =  0) 
and  at  the  model  top  the  vertical  velocity  has  a  zero  gradient. 
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3.  Model  results 


The  model  was  applied  to  four  cases,  i.e.,  a  single  building  in  the  center  of  the  grid  to  test  wind 
flow  symmetry,  an  array  of  five  buildings  (also  to  simulate  wind  flow  symmetry),  a  forest 
canopy  centered  in  the  model  grid  to  simulate  leading  edge  flow  separation,  and  a  larger  and 
more  complex  building  array  to  simulate  wind  flow  characteristics  around  the  IOL  and  adjacent 
ARL  structures.  Grid  spacing  for  the  single  and  multiple  building  array  was  unifonn,  i.e.,  Ax  = 
Ay  =  Az  =  1.0  m,  where  the  number  of  grid  points  was  27  x  27  x  12  for  the  single  building  and  50 
x  50  x  20  for  the  array.  Here,  At  =  0.001  s  and  the  model  was  executed  for  14K  time-steps.  For 
the  forest  application,  the  grid  spacing  was  Ax  =  Ay  =  2.0  m  and  Az  =  1.0  m,  where  the  number  of 
grid  points  was  60  x  60  x  30.  Here,  At  =  0.001  5  and  the  model  was  executed  for  40K  time-steps. 
For  the  IOL  and  adjacent  buildings,  the  grid  spacing  was  Ax  =  Ay  =  2.5  in  and  Az  =  2.0  in,  where 
the  number  of  grid  points  was  100  x  100  x  25.  Here,  in  contrast,  At  =  0.001  s  and  the  model  was 
executed  for  20K  time-steps.  The  computer  model  was  executed  on  a  1.8  GHz  desktop  PC  with 
approximately  512  MB  RAM.  Model  runtimes  were  approximately  2  minutes  for  the  single 
building,  20  minutes  for  the  five  building  array,  50  minutes  for  the  forest  canopy,  and  2  hours 
and  20  minutes  for  the  IOL  and  adjacent  buildings. 

Figure  1  presents  horizontal  and  vertical  cross-sections  of  the  vector  wind  field  around  a  single 
building  located  in  the  middle  of  the  model  grid.  Here,  the  initial  wind  velocities  are  u  = 

3.0  ms~l  and  v  =  w  =  0.0  ms~l.  The  horizontal  cross-section  is  taken  at  z  =  5.0  m  above  ground 
level.  This  result  is  very  encouraging.  The  wind  flow  patterns  going  around  (and  over  the  top 
of)  the  single  building  appear  quite  reasonable. 


Figure  1.  Horizontal  (a)  and  vertical  (b)  cross-section  of  the  vector  wind  field  around  a  single  building.  Initial  wind 
velocities  are  u  =  3.0  ms  and  v  =  w  =  0.0  ms1.  The  horizontal  cross-section  is  taken  at  z  =  5.0  m. 
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Figure  2  presents  horizontal  and  vertical  cross-sections  of  the  vector  wind  field  around  a  five 
building  array.  Again,  the  initial  wind  velocities  are  u  =  3.0  ms  1  and  v  =  w  =  0.0  ms  1  and  the 
cross-sections  are  taken  at y  =13.0  m  from  the  south  edge  andz  =  5.0  m  above  ground  level. 
These  results  are  also  quite  encouraging,  especially  since  they  provide  good  visual  confirmation 
(i.e.,  strong  symmetry)  that  the  program  coding  is  error-free. 


Figure  2.  Horizontal  (a)  and  vertical  (b)  cross-section  of  the  vector  wind  field  around  a  five  building  array.  The 
initial  wind  velocities  are  u  =  3.0  ms  and  v  =  w  =  0.0  ms1.  The  horizontal  cross-section  is  taken  at 
y  =  13.0  m  and  z  =  5.0  m. 


Figure  3  presents  vertical  cross-sections  of  the  wind  field  around  a  uniform  forest  canopy  located 
in  the  middle  of  the  model  grid.  The  canopy  height  is  h=  15.0  m.  The  leaf  area  index  is  LAI  =  4. 
In  addition,  a  generic  leaf  area  density  profile  is  assumed  (see  ref.  6,7).  Although  somewhat 
difficult  to  visualize  using  small  vector  arrows,  the  wind  field  (shown  in  figure  3a)  does  separate 
(diverge)  at  the  forest  edge.  Part  of  the  flow  goes  downward  into  the  trunk  space  of  the  forest, 
below  the  layer  of  leaves  and  branches,  and  part  of  the  flow  goes  upward  through  the  canopy  top 
and  above.  While  streamlines  may  be  generated  in  future  works  to  better  visualize  these  model 
data,  this  result  is  highly  consistent  with  previously  published  works  (16,  19,  and  20).  Moreover, 
figure  3b  shows  clearly  that  the  magnitude  of  the  winds  decrease  with  distance  into  the  forest 
(due  to  cumulative  drag  effects),  especially  in  the  upper  half  of  the  forest  stand.  In  addition, 
secondary,  low  level  wind  speed  maxima  develop  through  the  forest  leading  edge. 
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(b) 


Figure  3.  Vertical  cross-section  of  (a)  the  vector  wind  field  and  (b)  wind  velocity  contours  (units  ms  *)  around  a 
uniform  forest  canopy.  The  canopy  height  is  h  =  15.0  m.  The  leaf  area  index  is  LAI  =  4. 


Figure  4  presents  selected  wind  velocity  profiles  taken  from  the  model  data  shown  in  figure  3, 
i.e.,  the  initial  wind  velocity  (log-law)  profile,  the  profile  +20  m  upwind,  the  profile  at  the 
leading  edge,  the  profile  at  the  forest  center,  the  profile  at  the  trailing  edge,  and  the  profile  +20  m 
downwind  of  the  forest.  The  wind  speed  profile  in  the  open  fields  is  logarithmically  increasing 
with  height  above  the  roughness  plane.  Inside  the  forest,  the  winds  are  shown  to  decrease 
rapidly  as  momentum  becomes  depleted  through  the  layers  of  leaves  and  branches  (due  to 
increased  drag).  In  addition,  for  the  profile  in  the  forest  center,  the  model  produces  a  secondary 
wind  speed  maximum  at  a  height  of  about  z  =  4.0  m.  Shinn  ( 16)  and  Shaw  (21)  discuss  such 
low-level  wind  maxima  in  detail.  Note  that  the  model  results  shown  here  will  be  sensitive  to 
variations  in  the  assumed  leaf  area  index  (LAI)  and  leaf  area  density  distribution  in  the  canopy 
(6,7). 
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Figure  4.  Selected  wind  velocity  profdes  taken  from  the  results  shown  in  figure  3,  i.e.,  the  initial  wind  velocity  (log- 
law)  profile  (solid  line),  the  profile  +20  m  upwind  (large  dashed  line),  at  the  leading  edge  (dotted  line), 
at  the  forest  center  (dash-dotted  line),  at  the  trailing  edge  (dash-double  dotted  line),  and  +20  m 
downwind  of  the  forest  (small  dashed  line). 

Figure  5  presents  a  horizontal  cross-section  of  the  vector  wind  field  around  the  IOL  and  adjacent 
buildings.  Here,  the  initial  wind  velocities  are  u  =  w=  0.0  ms  1  and  v  =  4.0  ms~\  The  horizontal 
cross-section  is  taken  at  z  =  6.0  m  above  ground  level.  The  (dashed)  outlined  area  identifies 
forests  along  the  western  border  of  the  building  complex,  for  which  a  constant  drag  force  with 
height  (up  to  the  canopy  top,  h  =  18.0  m)  is  assumed.  More  complete  canopy  characterization 
data,  like  those  implemented  in  the  previous  example,  may  be  incorporated  in  future  works. 
Nevertheless,  this  model  result  is  very  appealing.  It  demonstrates  that  we  have  completed 
several  successful  steps  towards  obtaining  a  useful  mathematical  representation  of  the  wind  flow 
around  the  IOL  and  adjacent  buildings.  Naturally,  these  results  (as  well  as  those  shown  above) 
will  need  to  be  confirmed  by  measurements.  Yet,  if  determined  to  be  realistic,  this  kind  of  data 
will  be  very  important  for  electro-optic  and  acoustic  characterization  and  prediction. 
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Figure  5.  Horizontal  cross-section  of  the  vector  wind  field  around  the  IOL  and  adjacent  buildings.  The 

initial  wind  velocities  are  u  =  w  =  0.0  ms  1  and  v  =  4.0  ms1.  The  horizontal  cross-section  is  taken 
at  z  =  6.0  to  above  ground  level.  The  outlined  (dashed-double  dotted)  area  identifies  forests  along 
the  western  border  of  the  building  complex. 


4.  Discussion 


Note  that  the  wind  flow  model  presented  in  this  report  is  still  a  program  code  in  development, 
even  though  it  already  exhibits  some  powerful  capabilities.  Additional  algorithms  and/or 
alternate  numerical  schemes  may  be  implemented  in  future  works  to  improve  the  model.  For 
example,  temperature  and  moisture  prediction  may  be  investigated  via  the  following  expressions, 
i.e., 
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where  <9  is  air  temperature  (in  units  °K),  q  is  specific  humidity  (in  units  kg  kg'1),  Kg  and  Kq  are  the 
eddy  diffusivities  for  heat  and  water  vapor,  respectively  (in  units  ms'  ),  and  57;  and  Sq  are  the 
source  terms  for  heat  and  moisture.  The  quantity  Sg  for  forests  can  be  determined  via  the 
radiation  and  energy  budget  equations  presented  in  (7).  Alternately,  a  simpler  form  for  the  heat 
source  was  shown  by  Shaw  and  Schumann  (22),  i.e., 


dQ(z) 


dz 


(9) 


where  Q(z)  =  Q(h)exp(-aF).  Here,  Q  is  the  kinematic  heat  flux,  h  is  the  height  of  the  canopy 

h 

top,  a=  0.6  is  an  extinction  coefficient,  and  F  =  \  Adz  is  the  non-dimensional  cumulative  leaf 

Z 

area  index.  Naturally,  for  open  fields,  buildings,  and/or  pavement,  other  suitable  heat  source, 
moisture  source,  and  energy  budget  relations  can  be  applied. 


Overall,  we  expect  that  the  current  research  (when  completed)  will  be  quite  useful  to  improve 
optical  turbulence  intensity  ( Cn2 )  prediction  and  analysis  around  the  ARL  A  LOT  Facility.  Our 
work  will  also  contribute  towards  improved  estimates  of  other  key  signal  propagation 
parameters,  such  as  the  variance  in  angle-of-arrival  fluctuation  and  the  log-intensity  variance  of 
transmitted  electromagnetic  signals.  As  an  example,  for  propagation  across  distances  less  than 
or  equal  to  approximately  4  km,  Beland  (23)  gives  an  expression  for  the  variance  in  angle-of- 
arrival  fluctuation  for  spherical  waves  as, 


i 

2.914  D  3 


k2(z) 


dz . 


(10) 


for  D  »  (A  L )2 ,  where  D  (~  0.090  m,  possibly)  is  aperture  diameter,  A  (~  0.94  pm,  possibly)  is 
wavelength,  L  is  optical  path  length,  and  z  is  the  vertical  coordinate.  Analogous  expressions 
have  been  presented  for  the  log-intensity  (or  log-amplitude)  variance  of  transmitted 
electromagnetic  signals,  so  that  the  measure  of  the  path-averaged  Cn2  due  to  scintillation  (i.e., 
temporal  fluctuations)  for  a  spherical  wave  can  be  expressed  as, 


0.56 


k 6 


|C2(z) 


(L-z) 


dz 
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where  k  is  the  wave  propagation  constant  Thus,  if  improved  microphysical 

characterization  models  will  provide  better  estimates  of  Cn2  along  more  complex  optical  lines- 
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of-site,  then  better  estimates  of  beam  displacement,  (cr  2A  j,  and  intensity  fluctuations,  (a2x  j  ,  are 

the  kinds  improved  work  product  that  may  be  useful  in  sensor  performance  (and/or  sensor  data) 
analysis.  Note  that  for  model  evaluation  purposes,  the  ARL  ALOT  can  provide  time  series  of 
optical  turbulence  (scintillometer)  and  other  pertinent  data,  e.g.,  like  those  shown  in  figure  6. 


Figure  6.  Time  series  of  (a)  path-averaged  optical  turbulence  (scintillometer)  data  and  (b)  microphysical 
characterization  data  collected  at  the  ARL  A  LOT  Facility  on  12  December  2004. 

Finally,  a  useful  expression  for  the  refractive  index  structure  parameter  (24-26)  is, 


C2„  (z)  =  b  “TA 
£ 


f  dn^~ 


ydz  j 


(12) 


where  b  is  a  constant,  K/ ,  is  the  turbulent  exchange  coefficient  for  heat,  e  is  the  energy  dissipation 
rate,  and  dnldz  is  the  vertical  gradient  of  the  index  of  refraction  ( n ).  In  equation  12,  the  vertical 
gradient  of  the  index  of  refraction  for  visible  and  near-infrared  wavelengths  (0.36  to  3  jum)  can 
be  expressed  as  a  function  of  the  partial  derivatives  for  potential  temperature  ( dOldz )  and 
moisture  (specific  humidity)  ( dqtdz ),  i.e., 
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where  the  dispersion  fonnulas  (as  functions  of  wavelength  ',  i.e.,  a=  A  ')  are, 


M i( A)  =23. 7134+ 


6839.397  45.473 


130 -a2  38.9 -a2 


and 


M 2 (A )  =  64.873 1  +  0. 58058  a2  -0.007115  a4  +  0. 0008851  a6  ■ 


(14) 


(15) 


5.  Summary  and  Conclusions 


This  report  presented  a  new  finite-difference  computer  model  to  predict  the  surface  layer  wind 
flow  around  single  and  multiple  building  arrays  and  forest  canopies.  The  model  accounts  for 
advection,  the  pressure  gradient,  diffusion,  and  drag  forces  due  to  vegetation.  The  model  code  is 
computationally  efficient  and  extremely  flexible  with  regard  to  modifications  and  debugging. 
Our  initial  model  results  are  quite  encouraging.  Hence,  we  have  successfully  demonstrated  new 
capabilities  towards  obtaining  a  useful  mathematical  representation  of  microphysical 
characterization  data  around  the  IOL  and  adjacent  buildings.  In  conclusion,  we  expect  that  these 
kinds  of  simulation  codes  will  be  an  important  vehicle  for  investigating  optical  turbulence  intensity 
(< Cn2 )  and  related  laser-optics  propagation  effects  in  and  around  complex  enviromnents. 
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